/*
 * Sampling_ACM_FEM_Controller.cpp
 *
 *      Author: Jerome RD Soine and Christoph A Brand
 *      Institution: Schwarz goup, Institute for Theoretical Physics, Heidelberg University, Germany
 */


#include "Sampling_ACM_FEM_Controller.h"
#include "Sampling_ACM_FEM/StatisticsManager.h"
#include "ACM_FEM/OptimizationObject.h"
#include "include/core/worker/WorkerManager.h"
#include "include/core/worker/Worker.h"
#include "include/core/utility/VectorOperations.h"


Sampling_ACM_FEM_Controller::Sampling_ACM_FEM_Controller(string parameter_file, string DSF_file, string TVA_file, string VSF_file, unsigned int n_calcs)
{
	OptimizationObject* opti_object = new OptimizationObject(parameter_file);
	system(("mkdir " + opti_object->get_project_folder()).c_str());
	opti_object->save_parameters(opti_object->get_project_folder() + "/parameter_settings.log");

	//Init ForceTransmission shapes from file
	vector<Shape*> shapes = Ellipse::get_shapes_from_file(opti_object->get_shape_file());
	vector<Ellipse*> ellipses(shapes.size());
	for(unsigned int i=0; i<shapes.size(); i++) ellipses[i] = (Ellipse *) shapes[i];

	//Save experimental displacements in vtk format for visualization
//	Point<2> p1(opti_object->get_mesh_point1()(0),opti_object->get_mesh_point1()(1));
//	Point<2> p2(opti_object->get_mesh_point2()(0),opti_object->get_mesh_point2()(1));
//	FieldInterpolation2D * interpol = new FieldInterpolation2D(p1, p2, opti_object->get_degree_of_refinement(), opti_object->get_disp_file());
//	interpol->output(opti_object->get_project_folder() + "/displacements.vtk");
//
//	if(opti_object->is_preshifting_enabled())
//	{
//			//Shift mesh to reference (=nondisplaced) state
//			ActiveCableModel::createShiftedVTK(
//										opti_object->get_cell_mesh_file(),
//										opti_object->get_project_folder() + "/tmp.vtk",
//										interpol);
//			opti_object->set_cell_mesh_file(opti_object->get_project_folder() + "/tmp.vtk");
//
//			//Shift shapes to corresponding positions
//			for(unsigned int i=0; i<ellipses.size(); i++){
//				vector<double> tomove = ellipses[i]->get_center();
//				Vector<double> tmp (2);
//				Point<2> p(tomove[0], tomove[1]);
//				interpol->point_evaluation(p, tmp);
//				for(int j=0; j<10;j++){
//					tomove[0]=tomove[0]-tmp[0];
//					tomove[1]=tomove[1]-tmp[1];
//					Vector<double> backShift(2);
//					Point<2> controlP(tomove[0], tomove[1]);
//					interpol->point_evaluation(controlP, backShift);
//					controlP[0] += backShift[0];
//					controlP[1] += backShift[1];
//
//					Point<2> distance = controlP-p;
//					tmp[0] = distance[0];
//					tmp[1] = distance[1];
//				}
//				ellipses[i]->set_center(tomove);
//			}
//
//	}

	int dofnr = ActiveCableModel_FA::getFiberList(opti_object->get_cell_mesh_file()).size()+1;
	vector<double> solution (dofnr, 0);
	cout << "starting solution: \n";
	VectorOperations::cout_std_vector(solution);
	ActiveCableModel_FA * acm = new ActiveCableModel_FA(opti_object->get_cell_mesh_file(),opti_object->get_project_folder(),50.0,solution, ellipses);
	FEMCalc * fem = new FEMCalc(opti_object->get_young_mod(),opti_object->get_poisson_ratio(),opti_object->get_mesh_point1(),
			opti_object->get_mesh_point2(),opti_object->get_degree_of_refinement(),acm->getForceTransmissionShapes());
	L2Comparator * comp = new L2Comparator(new DataContainer(opti_object->get_disp_file()));

	unsigned int n_bins = 10;
	double network_tension = 0.112425;

	// Construct Statistics
	StatisticsManager* DSF_statistcs = new StatisticsManager(DSF_file,n_bins);
	StatisticsManager* TVA_statistcs = new StatisticsManager(TVA_file,n_bins);
	StatisticsManager* VSF_statistcs = new StatisticsManager(VSF_file,n_bins);

	// Construct Cell Model
	vector<int> fiber_type_list = acm->get_fiber_types();

	//
	unsigned int n_threads = opti_object->get_n_threads();
	WorkerManager::initialize(acm, fem, comp, n_threads);
	WorkerManager* workermanager = WorkerManager::getInstance();

	// Iterate over calculations
	vector<Worker*> worker_list;
	vector<double> parameters (acm->n_parameters(),network_tension);
	for(unsigned int i=0;i<n_threads;i++)
		worker_list.push_back(workermanager->getWorker(parameters,-1));
	for(unsigned int i=0;i<worker_list.size();i++)
		worker_list[i]->unlock();

	Worker* preWorker = workermanager->getWorker(parameters,-1);
	preWorker->start();
	preWorker->join();
	vector<double> init_state = preWorker->get_final_configuration();
	preWorker->unlock();

	unsigned int counter=0;
	ofstream output;
	output.open((opti_object->get_project_folder()+"/L2_stat.txt").c_str());
	output << "# L2 value list \n";
	output.close();
	do
	{
		if(counter<(n_calcs-n_threads))
		{
			for(unsigned int n=0;n<n_threads;n++){
				vector<double> params = generate_new_sample(fiber_type_list,DSF_statistcs,TVA_statistcs,VSF_statistcs,network_tension);
				worker_list[n] = workermanager->getWorker(params,init_state,-1);
				worker_list[n]->start();
				counter++;
			}
			for(unsigned int n=0;n<n_threads;n++){
				worker_list[n]->join();
				output.open((opti_object->get_project_folder()+"/L2_stat.txt").c_str(), ios::app);
				output << worker_list[n]->get_fitness() << "\n";
				output.close();
				worker_list[n]->unlock();
			}
		}
		else
		{
			for(unsigned int n=0;n<(n_calcs-n_threads);n++){
				vector<double> params = generate_new_sample(fiber_type_list,DSF_statistcs,TVA_statistcs,VSF_statistcs,network_tension);
				worker_list[n] = workermanager->getWorker(params,init_state,-1);
				worker_list[n]->start();
				counter++;
			}
			for(unsigned int n=0;n<(n_calcs-n_threads);n++){
				worker_list[n]->join();
				output.open((opti_object->get_project_folder()+"/L2_stat.txt").c_str(), ios::app);
				output << worker_list[n]->get_fitness() << "\n";
				output.close();
				worker_list[n]->unlock();
			}
		}
	}
	while(counter < n_calcs);
	cout << "finished" << endl;

	delete DSF_statistcs;
	delete TVA_statistcs;
	delete VSF_statistcs;

}

vector<double> Sampling_ACM_FEM_Controller::generate_new_sample(vector<int> fiber_type_list, StatisticsManager* dsf_stat, StatisticsManager* tva_stat, StatisticsManager* vsf_stat, double network_tension)
{
	//Initialize distribution
	vector<double> parameters (fiber_type_list.size()+1,network_tension);
	unsigned int n_dsf=0,n_vsf=0,n_tva=0;
	for(unsigned int fib=0;fib<fiber_type_list.size();fib++){
		if(fiber_type_list[fib]==ActiveCableModel::TYPE_DORSAL)
				{parameters[fib+1]=dsf_stat->pick_values_from_distribution(); n_dsf++;}
		else if(fiber_type_list[fib]==ActiveCableModel::TYPE_VENTRAL)
				{parameters[fib+1]=vsf_stat->pick_values_from_distribution(); n_vsf++;}
		else if (fiber_type_list[fib]==ActiveCableModel::TYPE_TVA)
				{parameters[fib+1]=tva_stat->pick_values_from_distribution(); n_tva++;}
		else if (fiber_type_list[fib]==ActiveCableModel::TYPE_EDGEBUNDLE)
				{parameters[fib+1]=vsf_stat->pick_values_from_distribution(); n_vsf++;}
	}
	cout << "N_SFs = " << n_dsf << "\t" << n_tva << "\t" << n_vsf << "\n";
	return parameters;
}
